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We present results from comprehensive event-driven (ED) simulations of nonlinear 
pattern formation in freely-evolving granular gases. In particular, we focus on the 
the morphologies of density and velocity fields in the inhomogeneous cooling state 
(ICS). We emphasize the strong analogy between the ICS morphologies and pattern 
formation in phase ordering systems with a globally conserved order parameter. 



There has been much recent interest in the properties of powders or granular materials 0, 
These materials are of great scientific and technological relevance and display properties 
intermediate to those of solids and fluids. Perhaps the most relevant feature of powders is 
that the grains dissipate energy on collision. Therefore, granular materials exhibit interesting 
dynamical properties only when the collisional energy loss is compensated by a continuous 
input of energy 

It is also relevant to investigate the dynamical evolution of a freely-evolving homogeneous 
system of inelastic granular particles. In the initial stages, the system continuously loses 
energy iri_^a homogeneous cooling state (HCS), where the density field is approximately 
uniform [3]. However, the HCS is unstable to density fluctuations, and the system evolves 
into an inhomogeneous cooling state (ICS), where particle-rich clusters are formed and 



grow 



There is a good understanding of the HCS and the instabilities which result 



in the ICS 
evolution of the ICS 



. However, there is only a limited understanding of the nonlinear 
3, 0|. In this letter, we study nonlinear domain growth processes 
for the granular density and velocity fields in the asymptotic time-regime. We present 
results from comprehensive event-driven (ED) simulations for a wide range of inelasticity 
and density parameters. In particular, we invoke analogies from studies of phase ordering 
dynamics P, to obtain a quantitative characterization of the evolving morphologies in 
the ICS. The similarities between freely-evolving granular gases and phase ordering systems 
lave been briefly discussed by Van Noije, Ernst and co-authors j^l, and Baldassari et al. 
^. Of course, we stress that the mechanisms driving pattern evolution are very different in 
both cases. 

We consider a system of identical inelastic hard spheres (with mass m = 1 and diameter 
a = 1) in a d = 2 box of size [NfyCry, with periodic boundary conditions. There is no 



external input of energy in our simulations. We use an event-driven algorithm |lll|, which 



only keeps track of particle collisions, to evolve the system. After a collision between the 
i*^ and j*^ particles, having velocities Vi and vj respectively, the new velocities are if^ = 
Vi — ^-^[n ■ {vi — Vj)]h, Vj = Vj + ■ {vi — Vj)]h. Here, h is the unit vector parallel to the 
relative position of the particles. For elastic collisions, we have the coefficient of restitution 
e = 1. For granular materials, e < 1 in general. 

The initial condition for each run consisted of particles with a homogeneous spatial 
distribution and a Maxwellian velocity distribution. For circular particles, the average pack- 
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ing fraction is = tiN / {ANl). We always fix A*";, = 256 and vary A^, as specified below. We 
will label our results using the number fraction, n = N/N^. Numerical data are obtained 
for the following parameter sets: (i) = 30000 (or n ~ 0.46, (p ~ 0.36), and e ranging from 
0.80 to 0.975 in steps of 0.025; (ii) e = 0.90, and ranging from 10000 to 40000 in steps 
of 5000, i.e., n ~ 0.15 to 0.61. We characterize the dynamical evolution of the granular gas 
using various statistical quantities, which are calculated as an average over 5 independent 
runs, i.e., correlation functions, structure factors, and domain growth laws for the density 
and velocity fields. 

For early times, inelastic collisions result in cooling as ^^^^ = — ^C(j(T)T, where T(t) is 
the temperature at time t; e = 1 — e^; and uj{T) is the collision frequency at temperature 
T j^. The form of uj(T) is approximately determined from the Enskog theory for elastic 
hard spheres as coiT) ~ co'(To)(T/To)"'^/^, where Tq is the granular temperature at t = 0. 

= TqC"^^, where 

r(t) is the average number of collisions per particle upto time t. Our subsequent results 
are presented in terms of r, as this constitutes a reasonable measure of time in the present 
context. The values of r are obtained directly from the ED simulations. We should stress 
that r(t) vs. t fiuctuates considerably for individual runs in the ICS regime. 

In Fig. 1, we show evolution snapshots for the d = 2 granular gas with n ~ 0.46. The LHS 
and RHS frames refer to the coarse-grained velocity field v{f,T) and density field ilj{r,T), 
respectively. The coarse-grained fields at a lattice point are obtained by averaging over boxes 
of size (5o")^ centered at that point. To clarify the nature of pattern formation, we have 
"hardened" the velocity field in Fig. 1, i.e., the length of all vectors has been set to unity. 
Points where the velocity field is zero, due to the absence of particles in the coarse-graining 
box, are unmarked in Fig. 1. The density field on the RHS of Fig. 1 is depicted in a binary 
representation. We introduce the order parameter %l){r,T) with values +1 (—1) at points 
where the number density is larger (less) than the average number density. 

The upper frames in Fig. 1 correspond to e = 0.85 and r = 150. The velocity field is 
characterized by the emergence and diffusive coarsening of vortices. There is a progressive 
parallelization of the local velocity field due to inelastic dissipation of the normal velocity 
components. However, as the total momentum is fixed at zero and the system is doubly- 
periodic, defects will be present in the velocity field for all times. Typically, pattern evolution 
in the density field is slower than that for the velocity field 0|. Nevertheless, well-defined 



This yields Haff's cooling law for the HCS [3| as T(t) = T[ 



1 I e^CTp) . 

2d ^ 



4 



clusters corresponding to the ICS are already seen at r = 50 (not shown here). These clusters 
grow with time in a manner reminiscent of phase ordering systems P]. The vortex centers 
in the velocity field are also shown (as black circles) in the frame on the RHS. The vortex 
density diminishes with time, and vortices are primarily confined to regions of low density, 
as there is a rapid parallelization of velocities in the high-density (solid-like) region due to 
multiple collision processes. There is no strong correlation between the location of defects 
in the velocity field (i.e., vortices) and defects in the density field (i.e., interfaces or domain 
boundaries). The middle frames in Fig. 1 correspond to e = 0.90. The broad features are 
the same as those for e = 0.85, though the time-scales are slower. All time-scales for pattern 
formation diverge as e ^ 1, i.e., the clustering instability is delayed for more elastic systems. 
The perfectly elastic case (e = 1) is singular, and does not exhibit clustering. 

A coarse-grained description of the granular fiuid is provided by nonlinear hydrodynamic 
equations for the density, velocity and temperature fields, in conjunction with an energy loss 
term Wakou et al. have shown that fiuctuations in these quantities o bey TDGL-like 
equations of phase ordering dynamics with a nonconserved order parameter 10] . Typically, 
these equations are obtained as the overdamped limit of a Hamiltonian formulation. In this 
letter, we clarify the TDGL-description of granular dynamics. 

The evolution of the density field is analogous to phase ordering dynamics in two- 
component (AB) mixtures, which is described by the equation jlo| : 

= (-V^)- [V^(f, t) - ^(f, tf + V2v^(f, t)] , (1) 

where %l){r,t) differentiates between A-rich {ip = +1) and i?-rich {tp = —1) regions. Eq. 
with m = is the TDGL equation, and describes nonconserved systems. In this case, the 
coarsening system is characterized by a diffusive growth law, L{t) ~ t^/^ jslQl, where L{t) 
is the characteristic domain size at time t. On the other hand, Eq. (Q) with m = 1 is the 
Cahn-Hilliard (CH) equation, and describes conserved systems. In this case, the domain 

nn 

growth process obeys the Lifshitz-Slyozov (LS) law, L(t) ~ t^^^ lallOl. Finally, there have 
also been studies of Eq. with m 0+, which is referred to as the globally-conserved (GC) 
TDGL equation jisl. This model is known to be in the same dynamical universality class 
as the nonconserved TDGL equation with the difference that the initial composition is 
preserved in the GC-TDGL equation. 

Similarly, the evolution of the velocity field in the ICS is comparable to ordering dynamics 
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in the XY model, which is the vector version of Eq. 



dip{f, t) 



j2\m 



tP{r,t)-\4j{f,t)\'i;{f,t)+V'i^{f,t) , (2) 



dt ' 

where ip = {'ipi,'4'2)- The evolution in Eq. (0) parallelizes ip locally via the annihilation of 
vortices and ant i- vortices, driven by a surface-tension reduction mechanism. Eq. (0) with 
m = corresponds to the case with nonconserved order parameter. The corresponding 
domain growth law 14j is Lv(t) ~ t^^^. On the other hand, Eq. (j2]) with m = 1 corresponds 
to the conserved XY model, where the dynamics locally conserves ip. 

The evolution of the granular fluid conserves both the density and velocity fields. Intu- 
itively, one may expect that the analogous phase ordering models would be the conserved 
(m = 1) versions of Eqs. (Q)-©- However, the evolution morphologies in the upper and 
middle frames of Fig. 1 are comparable to those for the case with nonconserved order pa- 
rameter . This is a consequence of the non-diffusive dynamics of granular particles, which 
move in straight lines until they collide with other particles. In this collision, the density 
and momentum are conserved quantities, of course. However, depending on the density of 
the fluid, the distance traveled by particles prior to collision may be considerable. There- 
fore, the variables are conserved on the macroscopic length-scale of the mean-free path, and 
not on the microscopic length-scale of the lattice spacing. Now, the ICS consists of regions 
of high density and low density. Granular particles stream relatively unhindered through 
the low-density regions and deposit on distant clusters. Typically, the conservation length 
scale is comparable to the length scale of the coarsening clusters, which diverges in time. 
Therefore, the density and velocity fields are globally conserved (GC), rather than locally 
conserved. To elucidate this analogy, the lower LHS frame in Fig. 1 shows the velocity field 
from the GC-XY model with (ip) = 0. The lower RHS frame in Fig. 1 shows the density 
field for the GC-TDGL model with {ip) = —0.08, corresponding to an average density of 
0.46. 

Before we quantify the morphologies in Fig. 1, we would like to study the HCS ICS 
crossover time. Brito and Ernst j3] use mode-coupling techniques to obtain the asymptotic 
energy decay (in the ICS) as 

X To (d-1 1 \ fAire \-^/^ 
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where ^± ~ \/~^o and ^ii ~ —Iq, with Iq being the time-independent mean-free path. A 



II 

comparison of Haff 's law and Eq. yields the crossover time Tc as the solution of 

r5V^^^^^(g)'xHV-\ (4) 

where cr = 1, and Vld = 2tc'^^'^ /T(d/2) is the d-dimensional solid angle. We consider the case 
e ^ 0, so that ^|| ^ For d = 2, Eq. (0]) simplifies as TcC'^'^" ~ nx(n)^, where we ignore 
prefactors. We use the Verlet-Levesque approximation for the d = 2 hard-sphere correlation 
function at contact, viz., = (1 — 70/16)(l — 0)^^. In Fig. 2(a), we plot e^'^" vs. Tc 

for n ~ 0.46 and a range of e-values. In Fig. 2(b), we plot Tce"^"^" vs. nx{nY for e = 0.90 
and a range of n- values. These approximately linear plots confirm the validity of the scaling 
behavior of Tc{e,n), as expected from Eq. 

Next, we examine the time- dependent structure factors of the density and velocity fields 
in the ICS. The evolving morphologies are characterized by unique length scales, and we 
expect the structure factors to exhibit dynamical scaling, i.e., S(k, t) = lijYfikt), where 
i is the relevant length scale and /(x) is a scaling function [lO|. In Fig. 3(a), we plot the 
scaled structure factors of the density field for (n, e) = (0.46,0.85), i.e., \n.[S^^{k,T){k)'^] 
vs. ln{k/{k)) from three different times (r ^ Tc). The quantity (k) is the first moment 
of the structure factor, and is related to the length scale as (k)^^ ~ L{t). The reasonable 
data collapse confirms the validity of dynamical scaling. The dashed line in Fi^^3(a) 
denotes the Fourier transform of the analytic result due to Ohta et al. (OJK) [16] for 
the correlation function of the nonconserved TDGL equation: g[x) = ^sin^^(7), where 
7 = e~^^,x = r/L. The OJK function provides an excellent description of our numerical 
data. The dot-dashed line in Fig. 3(a) is obtained from a numerical simulation of the CH 
equation with (ip) = —0.08. The CH result is at variance with our numerical results except 
for the tail region. In particular, the local conservation law dictates that Sc}i{k,t) ~ k'^ as 
However, the ICS structure factor appears to decay almost monotonically from 
A; = 0"^, as is usual for structure factors in ordering problems characterized by a nonconserved 
order parameter. (Of course, global density conservation dictates that S^^{0,t) = 0.) Fig. 
3(b) is analogous to 3(a), but corresponds to the case with e = 0.90. 

We have confirmed numerically that the scaled structure factor for the density field in 
the ICS has only a weak dependence on the off-criticality which measures the average 
density. Furthermore, it is numerically comparable to the OJK function for a broad range 
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of (z/')-values. The nature of defects in the ordering field dictates some general properties of 
the structure factor. For example, scattering off interfaces in the density field gives rise to a 
power-law (or Porod) decay in the structure factor tail, S^^{k,T) ~ f^j. large k 

which is seen clearly in Figs. 3(a)-(b). 

In Figs. 4(a)-(b), we plot scaled structure factors, \n.[Sm{k^T){k)'^] vs. ln(fc/(A;)), for the 
velocity field in the ICS. The dashed line denotes the Fourier transform of the Bray-Puri- 
Toyoki (BPT) function for ordering in the nonconserved XY model 



, , , n-i \^fn + l IM^ „ /I 1 n + 2 



^(i,i;^;7^), (5) 



where 7 = , x = r/L^; B{x, y) is the beta function; and F{a, b; c; z) is the hypergeometric 
function. The BPT result corresponds to defects with 0{n) symmetry, and the case with 
n = 2 is relevant here. Again, we see that our numerical data is described well by the BPT 
function, and is completely different from the result for the conserved XY model (obtained 
numerically, and denoted as a dot-dashed line in Figs. 4(a)-(b)). 

We have seen that vortex defects characterize the morphology of the ICS velocity field. 
These vortices give rise to a power-law decay for the structure factor tail, Syy{k, r) ~ j^-{d+2) 
for large k. The general BPT result in Eq. (jSJ exhibits a power-law tail, S{k,t) ~ k^^^^^\ 
and is relevant for the ICS in granular systems with d > 2. 

For the time-regimes considered here, a description in terms of two uncoupled order 
parameters is reasonable. A more complete description of the asymptotic ordering dynamics 
should account for holes in the velocity field due to absence of particles in some regions of 
space. This is provided by a model with spin- vacancy phase separation in conjunction with 
XY-like ordering of the spin variable. In the current context, the relevant model has coupled 
nonconserved (or globally conserved) dynamics for the two ordering fields. At a later stage, 
we will discuss in detail the applicability of this model to the ICS. 

Finally, we would like to briefly discuss domain growth laws in the ICS. We should stress 
that late-stage dynamics is affected by inelastic collapse, where a group of particles undergo 
an infinite number of collisions in a finite time ^|. The results discussed here correspond 
to the regime prior to inelastic collapse. Subsequent to the HCS ICS crossover, our 
numerical results (not shown here) are consistent with diffusive growth, L{t)^L^{t) ~ r^/^ 
0]. Furthermore, our numerical data is consistent with an asymptotic power-law behavior, 
r(t) ~ t^/'^, for a wide range of parameter values. Therefore, in real time, the length scales 
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behave as L(t), Li,{t) ~ t^^^ [t^. The domain growth law for the ICS density field is the same 
as the LS growth law for phase separation of binary mixtures. 

In conclusion, we have undertaken comprehensive ED simulations of nonlinear pattern 
formation in the density and velocity fields of an inelastic granular gas. We find that 
there is a close analogy between ICS morphologies and phase ordering systems with 
globally conserved order parameters. The nature of defects in the ordering fields dictates 
general properties of the relevant structure factors and correlation functions. Furthermore, 
our numerical results for the structure factors of the ICS density and velocity fields are 
described well by analytic results for nonconserved phase ordering systems. Clearly, the 
general formalism of phase ordering dynamics is of great utility in diverse problems of 
pattern formation as many features of evolution morphologies are determined by general 
principles, e.g., defect structures, conservation laws, etc. 
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Figure Captions 

Figure 1: Evolution pictures for the coarse-grained density and velocity fields in the ICS. 
The upper RHS frame shows the density field at r = 150 for (n, e) = (0.46, 0.85). Regions 
with density greater than the average are marked in grey. The black circles denote the 
vortex centers of the velocity field. The upper LHS frame shows the direction vectors of 
the corresponding velocity field. For clarity, we only show a 32^ corner (as depicted in the 
RHS frame) of the 256^ lattice. The middle frames are the density and velocity fields at 
T = 150 for (n, e) = (0.46,0.90). The lower RHS frame shows an evolution picture (at 
t — 25) for the GC-TDGL equation with (■0) = —0.08, corresponding to average density 
0.46. The lower LHS frame shows an evolution picture (at t — 25) for the GC-XY model 
with (tf) = 0. 

Figure 2: Dependence of Tc on system parameters. The crossover time is obtained 
from energy-decay plots (not shown here) as the point of deviation from Haff's law. The 
error bars are defined by the symbol sizes, (a) Plot of ei'^" vs. Tc for n ~ 0.46; and e ranging 
from 0.80 to 0.95. The solid line denotes the best linear fit to the data, (b) Plot of TcC'^'^" 
vs. nx{nY for e = 0.90; and n ranging from 0.08 to 0.53. 

Figure 3: (a) Scaled structure factors for the coarse-grained density field for 
(n,e) = (0.46,0.85). Data is plotted for t(> Tc) = 100,150,250. (b) Analogous to 
(a), but for (n, e) = (0.46,0.90). 

Figure 4: Analogous to Fig. 3, but for the coarse-grained velocity field. 
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